Code
library(here)
library(flextable)
library(plotly)
source(here("02_scripts/03a_fun_calib_td.R"))
source(here("02_scripts/06_interventions_models.R"))
theme_set(theme_minimal())CBA results point estimates
Scenario 1: Accounts for the effect of Fentanyl contamination starting from 2018 Scenario 2: Scenario 2: Doesn’t account for Fentanyl contamination and has the same probability from Illicit use to OD and from OD to OD death
Code
trace_tbl <- as_tibble(mod_basecase$m_M) %>%
mutate(cycle_num = 0:n_cycles) %>%
pivot_longer(1:n_states, names_to = "state", values_to = "count") %>%
mutate(mod = "No Interventions") %>%
bind_rows(., as_tibble(mod_nalx_1$m_M) %>%
mutate(cycle_num = 0:n_cycles) %>%
pivot_longer(1:n_states, names_to = "state", values_to = "count") %>%
mutate(mod = "Naloxone")) %>%
bind_rows(., as_tibble(mod_ss_1$m_M) %>%
mutate(cycle_num = 0:n_cycles) %>%
pivot_longer(1:n_states, names_to = "state", values_to = "count") %>%
mutate(mod = "Safer Supply")) %>%
bind_rows(., as_tibble(mod_pres_guid_1$m_M) %>%
mutate(cycle_num = 0:n_cycles) %>%
pivot_longer(1:n_states, names_to = "state", values_to = "count") %>%
mutate(mod = "Prescription Guidelines")) %>%
bind_rows(., as_tibble(mod_all_int_1$m_M) %>%
mutate(cycle_num = 0:n_cycles) %>%
pivot_longer(1:n_states, names_to = "state", values_to = "count") %>%
mutate(mod = "All Interventions"))
trace_tbl$mod <- factor(trace_tbl$mod,
levels = c("No Interventions", "Naloxone", "Safer Supply", "Prescription Guidelines", "All Interventions"),
labels = c("No Interventions", "Naloxone", "Safer Supply", "Prescription Guidelines", "All Interventions"))
v_extra_sevbi <- c(mod_basecase$extra_sev_bi, mod_nalx_1$extra_sev_bi,
mod_ss_1$extra_sev_bi, mod_pres_guid_1$extra_sev_bi,
mod_all_int_1$extra_sev_bi)
v_extra_modbi <- c(mod_basecase$extra_mod_bi, mod_nalx_1$extra_mod_bi,
mod_ss_1$extra_mod_bi, mod_pres_guid_1$extra_mod_bi,
mod_all_int_1$extra_mod_bi)
v_extra_od <- c(mod_basecase$extra_od_deaths, mod_nalx_1$extra_od_deaths,
mod_ss_1$extra_od_deaths, mod_pres_guid_1$extra_od_deaths,
mod_all_int_1$extra_od_deaths)
mod_names <- c("No Interventions", "Naloxone", "Safer Supply", "Prescription Guidelines", "All Interventions")
mod_names <- factor(mod_names, levels = mod_names, labels = mod_names)
for (j in 1:length(mod_names)){
trace_tbl$count[trace_tbl$state == "BO_OD_DEATH" &
trace_tbl$cycle_num == 180 &
trace_tbl$mod == mod_names[j]] <- trace_tbl$count[trace_tbl$state == "BO_OD_DEATH" &
trace_tbl$cycle_num == 180 &
trace_tbl$mod == mod_names[j]] + v_extra_od[j]
trace_tbl$count[trace_tbl$state == "BO_MOD_BI" &
trace_tbl$cycle_num == 180 &
trace_tbl$mod == mod_names[j]] <- trace_tbl$count[trace_tbl$state == "BO_MOD_BI" &
trace_tbl$cycle_num == 180 &
trace_tbl$mod == mod_names[j]] + v_extra_modbi[j]
trace_tbl$count[trace_tbl$state == "BO_SEVERE_BI" &
trace_tbl$cycle_num == 180 &
trace_tbl$mod == mod_names[j]] <- trace_tbl$count[trace_tbl$state == "BO_SEVERE_BI" &
trace_tbl$cycle_num == 180 &
trace_tbl$mod == mod_names[j]] + v_extra_sevbi[j]
}
trace_tbl$state <- factor(trace_tbl$state,
levels = v_state_names,
labels = v_state_names)Code
trace_tbl_2 <- as_tibble(mod_basecase_2$m_M) %>%
mutate(cycle_num = 0:n_cycles) %>%
pivot_longer(1:n_states, names_to = "state", values_to = "count") %>%
mutate(mod = "No Interventions") %>%
bind_rows(., as_tibble(mod_nalx_2$m_M) %>%
mutate(cycle_num = 0:n_cycles) %>%
pivot_longer(1:n_states, names_to = "state", values_to = "count") %>%
mutate(mod = "Naloxone")) %>%
bind_rows(., as_tibble(mod_ss_2$m_M) %>%
mutate(cycle_num = 0:n_cycles) %>%
pivot_longer(1:n_states, names_to = "state", values_to = "count") %>%
mutate(mod = "Safer Supply")) %>%
bind_rows(., as_tibble(mod_pres_guid_2$m_M) %>%
mutate(cycle_num = 0:n_cycles) %>%
pivot_longer(1:n_states, names_to = "state", values_to = "count") %>%
mutate(mod = "Prescription Guidelines")) %>%
bind_rows(., as_tibble(mod_all_int_2$m_M) %>%
mutate(cycle_num = 0:n_cycles) %>%
pivot_longer(1:n_states, names_to = "state", values_to = "count") %>%
mutate(mod = "All Interventions"))
trace_tbl_2$mod <- factor(trace_tbl_2$mod,
levels = c("No Interventions", "Naloxone", "Safer Supply", "Prescription Guidelines", "All Interventions"),
labels = c("No Interventions", "Naloxone", "Safer Supply", "Prescription Guidelines", "All Interventions"))
v_extra_sevbi_2 <- c(mod_basecase_2$extra_sev_bi, mod_nalx_2$extra_sev_bi,
mod_ss_2$extra_sev_bi, mod_pres_guid_2$extra_sev_bi,
mod_all_int_2$extra_sev_bi)
v_extra_modbi_2 <- c(mod_basecase_2$extra_mod_bi, mod_nalx_2$extra_mod_bi,
mod_ss_2$extra_mod_bi, mod_pres_guid_2$extra_mod_bi,
mod_all_int_2$extra_mod_bi)
v_extra_od_2 <- c(mod_basecase_2$extra_od_deaths, mod_nalx_2$extra_od_deaths,
mod_ss_2$extra_od_deaths, mod_pres_guid_2$extra_od_deaths,
mod_all_int_2$extra_od_deaths)
for (j in 1:length(mod_names)){
trace_tbl_2$count[trace_tbl_2$state == "BO_OD_DEATH" &
trace_tbl_2$cycle_num == 180 &
trace_tbl_2$mod == mod_names[j]] <- trace_tbl_2$count[trace_tbl_2$state == "BO_OD_DEATH" &
trace_tbl_2$cycle_num == 180 &
trace_tbl_2$mod == mod_names[j]] + v_extra_od_2[j]
trace_tbl_2$count[trace_tbl_2$state == "BO_MOD_BI" &
trace_tbl_2$cycle_num == 180 &
trace_tbl_2$mod == mod_names[j]] <- trace_tbl_2$count[trace_tbl_2$state == "BO_MOD_BI" &
trace_tbl_2$cycle_num == 180 &
trace_tbl_2$mod == mod_names[j]] + v_extra_modbi_2[j]
trace_tbl_2$count[trace_tbl_2$state == "BO_SEVERE_BI" &
trace_tbl_2$cycle_num == 180 &
trace_tbl_2$mod == mod_names[j]] <- trace_tbl_2$count[trace_tbl_2$state == "BO_SEVERE_BI" &
trace_tbl_2$cycle_num == 180 &
trace_tbl_2$mod == mod_names[j]] + v_extra_sevbi_2[j]
}
trace_tbl_2$state <- factor(trace_tbl_2$state,
levels = v_state_names,
labels = v_state_names)
trace_tbl <- trace_tbl %>%
mutate(scenario = "Scenario 1") %>%
bind_rows(., trace_tbl_2 %>%
mutate(scenario = "Scenario 2"))1 Cumulative OD Deaths
1.0.0.1 Cumulative OD Deaths - Scenario 1
Code
p1 <- ggplot(trace_tbl %>%
filter(state == "BO_OD_DEATH") %>%
filter(cycle_num != 0) %>%
filter((cycle_num + 1) %% 12 == 1) %>%
filter(scenario == "Scenario 1") %>%
mutate(year = rep(c(2015:2029), length(mod_names))),
aes(x = year, y = count, colour = mod)) +
geom_line() +
geom_point(size = 1) +
ylab("Cumulative OD Deaths") + xlab("Year") +
labs(colour = "Intervention")
ggplotly(p1)1.0.0.2 Cumulative OD Deaths - Scenario 2
Code
p2 <- ggplot(trace_tbl %>%
filter(state == "BO_OD_DEATH") %>%
filter(cycle_num != 0) %>%
filter((cycle_num + 1) %% 12 == 1) %>%
filter(scenario == "Scenario 2") %>%
mutate(year = rep(c(2015:2029), length(mod_names))),
aes(x = year, y = count, colour = mod)) +
geom_line() +
geom_point(size = 1) +
ylab("Cumulative OD Deaths") + xlab("Year") +
labs(colour = "Intervention")
ggplotly(p2)1.0.0.3 Both scenarios
Code
plotly::ggplotly(ggplot(trace_tbl %>%
filter(state == "BO_OD_DEATH") %>%
filter(cycle_num != 0) %>%
filter((cycle_num + 1) %% 12 == 1) %>%
mutate(year = rep(c(2015:2029), length(mod_names)*2)),
aes(x = year, y = count, colour = mod)) +
geom_line() +
geom_point(size = 1) +
ylab("Cumulative OD Deaths") + xlab("Year") +
labs(colour = "Intervention") +
facet_wrap(~scenario) +
labs(caption = "Scenario 1: Accounts for the effect of Fentanyl contamination starting from 2018 \n
Scenario 2: Doesn't account for Fentanyl contamination and has the same probability from \n Illicit use to OD and from OD to OD death"))2 New OD Deaths per year
Code
inci_od_deaths_basecase <- inci_od_deaths_nalox <- inci_od_deaths_ss <- inci_od_deaths_pg <- inci_od_deaths_all <- inci_od_deaths_basecase_2 <- inci_od_deaths_nalox_2 <- inci_od_deaths_ss_2 <- inci_od_deaths_pg_2 <- inci_od_deaths_all_2 <- rep(NA, 14)
for (i in 1:14){
yr <- 2015 + i
inci_od_deaths_basecase[i] <- mod_basecase$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_OD_DEATH"] -
mod_basecase$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr - 1 &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_OD_DEATH"]
inci_od_deaths_basecase_2[i] <- mod_basecase_2$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_OD_DEATH"] -
mod_basecase_2$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr - 1 &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_OD_DEATH"]
inci_od_deaths_nalox[i] <- mod_nalx_1$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_OD_DEATH"] -
mod_nalx_1$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr - 1 &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_OD_DEATH"]
inci_od_deaths_nalox_2[i] <- mod_nalx_2$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_OD_DEATH"] -
mod_nalx_2$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr - 1 &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_OD_DEATH"]
inci_od_deaths_ss[i] <- mod_ss_1$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_OD_DEATH"] -
mod_ss_1$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr - 1 &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_OD_DEATH"]
inci_od_deaths_ss_2[i] <- mod_ss_2$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_OD_DEATH"] -
mod_ss_2$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr - 1 &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_OD_DEATH"]
inci_od_deaths_pg[i] <- mod_pres_guid_1$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_OD_DEATH"] -
mod_pres_guid_1$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr - 1 &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_OD_DEATH"]
inci_od_deaths_pg_2[i] <- mod_pres_guid_2$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_OD_DEATH"] -
mod_pres_guid_2$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr - 1 &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_OD_DEATH"]
inci_od_deaths_all[i] <- mod_all_int_1$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_OD_DEATH"] -
mod_all_int_1$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr - 1 &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_OD_DEATH"]
inci_od_deaths_all_2[i] <- mod_all_int_2$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_OD_DEATH"] -
mod_all_int_2$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr - 1 &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_OD_DEATH"]
}
inc_od_death_tbl <- data.frame(intervention = c(rep(mod_names, each = 14)),
year = c(rep(2016:2029, 5)),
inci_od_deaths = c(inci_od_deaths_basecase,
inci_od_deaths_nalox,
inci_od_deaths_ss,
inci_od_deaths_pg,
inci_od_deaths_all),
scenario = rep("Scenario 1", 14*5)) %>%
bind_rows(., data.frame(intervention = c(rep(mod_names, each = 14)),
year = c(rep(2016:2029, 5)),
inci_od_deaths = c(inci_od_deaths_basecase_2,
inci_od_deaths_nalox_2,
inci_od_deaths_ss_2,
inci_od_deaths_pg_2,
inci_od_deaths_all_2),
scenario = rep("Scenario 2", 14*5)))
saveRDS(inc_od_death_tbl, file = here("01_data/inc_od_death_tbl_mod.RDS"))
target_od_deaths <- calib_target_tbl %>%
filter(group == "total_od_deaths") %>%
select(year, target) %>%
rename(inci_od_deaths = target) %>%
mutate(intervention = "Target")
saveRDS(target_od_deaths, file = here("01_data/inc_od_death_tbl_target.RDS"))
inc_od_death_tbl$intervention <- factor(inc_od_death_tbl$intervention,
levels = mod_names, labels = mod_names)2.0.0.1 Scenario 1
Code
p3 <- ggplot(inc_od_death_tbl %>% filter(scenario == "Scenario 1"),
aes(x = year, y = inci_od_deaths, color = intervention)) +
geom_line() +
geom_point(size = 1) +
geom_point(data = target_od_deaths, aes (x = year, y = inci_od_deaths)) +
ylab("New OD Deaths per year") + xlab("Year") +
labs(colour = "Intervention") +
scale_color_brewer(palette = "Dark2")
ggplotly(p3)2.0.0.2 Scenario 2
Code
p4 <- ggplot(inc_od_death_tbl %>% filter(scenario == "Scenario 2"),
aes(x = year, y = inci_od_deaths, color = intervention)) +
geom_line() +
geom_point(size = 1) +
geom_point(data = target_od_deaths, aes (x = year, y = inci_od_deaths)) +
ylab("New OD Deaths per year") + xlab("Year") +
labs(colour = "Intervention") +
scale_color_brewer(palette = "Dark2")
ggplotly(p4)2.0.0.3 Both scenarios
Code
plotly::ggplotly(ggplot(inc_od_death_tbl,
aes(x = year, y = inci_od_deaths, color = intervention)) +
geom_line() +
geom_point(size = 1) +
geom_point(data = target_od_deaths, aes (x = year, y = inci_od_deaths)) +
ylab("New OD Deaths per year") + xlab("Year") +
labs(colour = "Intervention") +
scale_color_brewer(palette = "Dark2") +
facet_wrap(~scenario))3 Cumulative Deaths
3.0.0.1 Scenario 1
Code
p5 <- ggplot(trace_tbl %>%
filter(state == "BO_DEATH") %>%
filter(cycle_num != 0) %>%
filter((cycle_num + 1) %% 12 == 1) %>%
filter(scenario == "Scenario 1") %>%
mutate(year = rep(c(2015:2029), length(mod_names))),
aes(x = year, y = count, colour = mod)) +
geom_line() +
ylab("Cumulative Deaths (excluding OD-deaths)") + xlab("Year") +
labs(colour = "Intervention")
ggplotly(p5)3.0.0.2 Scenario 2
Code
p6 <- ggplot(trace_tbl %>%
filter(state == "BO_DEATH") %>%
filter(cycle_num != 0) %>%
filter((cycle_num + 1) %% 12 == 1) %>%
filter(scenario == "Scenario 2") %>%
mutate(year = rep(c(2015:2029), length(mod_names))),
aes(x = year, y = count, colour = mod)) +
geom_line() +
ylab("Cumulative Deaths (excluding OD-deaths)") + xlab("Year") +
labs(colour = "Intervention")
ggplotly(p6)3.0.0.3 Both scenarios
Code
plotly::ggplotly(ggplot(trace_tbl %>%
filter(state == "BO_DEATH") %>%
filter(cycle_num != 0) %>%
filter((cycle_num + 1) %% 12 == 1) %>%
mutate(year = rep(c(2015:2029), length(mod_names)*2)),
aes(x = year, y = count, colour = mod)) +
geom_line() +
ylab("Cumulative Deaths (excluding OD-deaths)") + xlab("Year") +
labs(colour = "Intervention") +
facet_wrap(~scenario))4 Trace Plots
Code
set.seed(122345)
colors_pal <- Polychrome::createPalette(n_states,
c("#084C61", "#DB504A",
"#E3B505", "#4F6D7A",
"#56A3A6"))
names(colors_pal) <- NULL
ggplotly(ggplot(trace_tbl %>%
filter(grepl("BPO", state)), aes(x = cycle_num, y = count, colour = state)) +
geom_line() +
ylab("Count") + xlab("Cycle Number") +
ggrepel::geom_text_repel(data = trace_tbl %>% filter(cycle_num == 180) %>%
filter(grepl("BPO", state)),
aes(x = cycle_num, y = count, label = state),
check_overlap = T, size = 3) +
scale_color_manual(values = colors_pal) +
# facet_wrap(~mod) +
facet_grid(rows = vars(mod), cols = vars(scenario),
switch = "y"))Code
plotly::ggplotly(ggplot(trace_tbl %>%
filter(state %in% c("BPO_MISUSE", "BI_ILLICIT",
"BS_OAT_INI", "BS_OAT_MAINT")), aes(x = cycle_num, y = count, colour = state)) +
geom_line() +
ylab("Count") + xlab("Cycle Number") +
ggrepel::geom_text_repel(data = trace_tbl %>%
filter(state %in% c("BPO_MISUSE", "BI_ILLICIT",
"BS_OAT_INI", "BS_OAT_MAINT")) %>% filter(cycle_num == 180),
aes(x = cycle_num, y = count, label = state),
check_overlap = T, size = 3) +
scale_color_manual(values = colors_pal) +
# facet_wrap(~mod) +
facet_grid(rows = vars(mod), cols = vars(scenario),
switch = "y"))Code
plotly::ggplotly(ggplot(trace_tbl %>%
filter(state %in% c("BS_OAT_SS", "BS_SS")), aes(x = cycle_num, y = count, colour = state)) +
geom_line() +
ylab("Count") + xlab("Cycle Number") +
ggrepel::geom_text_repel(data = trace_tbl %>%
filter(state %in% c("BS_OAT_SS", "BS_SS")) %>% filter(cycle_num == 180),
aes(x = cycle_num, y = count, label = state),
check_overlap = T, size = 3) +
scale_color_manual(values = colors_pal) +
# facet_wrap(~mod) +
facet_grid(rows = vars(mod), cols = vars(scenario),
switch = "y"))Code
plotly::ggplotly(ggplot(trace_tbl %>%
filter(state %in% c("BO_OD_ILLICIT", "BO_OD_RX",
"BO_MOD_BI", "BO_SEVERE_BI")), aes(x = cycle_num, y = count, colour = state)) +
geom_line() +
ylab("Count") + xlab("Cycle Number") +
ggrepel::geom_text_repel(data = trace_tbl %>%
filter(state %in% c("BO_OD_ILLICIT", "BO_OD_RX",
"BO_MOD_BI", "BO_SEVERE_BI")) %>% filter(cycle_num == 180),
aes(x = cycle_num, y = count, label = state),
check_overlap = T, size = 3) +
scale_color_manual(values = colors_pal) +
# facet_wrap(~mod) +
facet_grid(rows = vars(mod), cols = vars(scenario),
switch = "y"))Code
plotly::ggplotly(ggplot(trace_tbl %>%
filter(state %in% c("BR_ILLICIT", "BR_OAT_INI",
"BR_OAT_MAINT", "BR_OAT_SS",
"BR_SS", "BR_OD_ILLICIT")), aes(x = cycle_num, y = count, colour = state)) +
geom_line() +
ylab("Count") + xlab("Cycle Number") +
ggrepel::geom_text_repel(data = trace_tbl %>%
filter(state %in% c("BR_ILLICIT", "BR_OAT_INI",
"BR_OAT_MAINT", "BR_OAT_SS",
"BR_SS", "BR_OD_ILLICIT")) %>% filter(cycle_num == 180),
aes(x = cycle_num, y = count, label = state),
check_overlap = T, size = 3) +
scale_color_manual(values = colors_pal) +
# facet_wrap(~mod) +
facet_grid(rows = vars(mod), cols = vars(scenario),
switch = "y"))5 CBA results
5.0.1 Scenario 1
Code
cost_diff <- c(mod_nalx_1$total_net_present_cost - mod_basecase$total_net_present_cost,
mod_ss_1$total_net_present_cost - mod_basecase$total_net_present_cost,
mod_pres_guid_1$total_net_present_cost - mod_basecase$total_net_present_cost,
mod_all_int_1$total_net_present_cost - mod_basecase$total_net_present_cost)
cost_diff_per <- c(round((cost_diff/c(mod_basecase$total_net_present_cost))*100,2))
deaths_eff_tbl <- trace_tbl %>%
filter(scenario == "Scenario 1") %>%
filter(cycle_num == 180) %>%
select(-cycle_num) %>%
pivot_wider(names_from = mod, values_from = count) %>%
filter(state %in% c("BO_DEATH", "BO_OD_DEATH")) %>%
mutate(nalx_effect = Naloxone - `No Interventions`,
ss_effect = `Safer Supply` - `No Interventions`,
pg_effect = `Prescription Guidelines` - `No Interventions`,
all_effect = `All Interventions` - `No Interventions`,
nalx_effect_per = round(((Naloxone - `No Interventions`)/`No Interventions`)*100, 2),
ss_effect_per = round(((`Safer Supply` - `No Interventions`)/`No Interventions`)*100, 2),
pg_effect_per = round(((`Prescription Guidelines` - `No Interventions`)/`No Interventions`)*100, 2),
all_effect_per = round(((`All Interventions` - `No Interventions`)/`No Interventions`)*100, 2)) %>%
select(-c(Naloxone, `No Interventions`, `Safer Supply`, `Prescription Guidelines`, `All Interventions`))
total_death_oddeaths <- trace_tbl %>%
filter(scenario == "Scenario 1") %>%
filter(cycle_num == 180) %>%
select(-cycle_num) %>%
pivot_wider(names_from = mod, values_from = count) %>%
filter(state %in% c("BO_DEATH", "BO_OD_DEATH"))
od_deaths_diff <- c(deaths_eff_tbl$nalx_effect[1], deaths_eff_tbl$ss_effect[1],
deaths_eff_tbl$pg_effect[1], deaths_eff_tbl$all_effect[1])
od_deaths_diff_per <- c(deaths_eff_tbl$nalx_effect_per[1], deaths_eff_tbl$ss_effect_per[1],
deaths_eff_tbl$pg_effect_per[1], deaths_eff_tbl$all_effect_per[1])
deaths_diff <- c(deaths_eff_tbl$nalx_effect[2], deaths_eff_tbl$ss_effect[2],
deaths_eff_tbl$pg_effect[2], deaths_eff_tbl$all_effect[2])
deaths_diff_per <- c(deaths_eff_tbl$nalx_effect_per[2], deaths_eff_tbl$ss_effect_per[2],
deaths_eff_tbl$pg_effect_per[2], deaths_eff_tbl$all_effect_per[2])
saveRDS(cbind.data.frame(cost_diff, cost_diff_per,
od_deaths_diff, od_deaths_diff_per,
deaths_diff, deaths_diff_per), here("01_data/point_estiamte.RDS"))
saveRDS(total_death_oddeaths, here("01_data/total_deaths_oddeaths.RDS"))
costs <- paste0(round(cost_diff/1000000, 0), " (", cost_diff_per, "%)")
deaths <- paste0(round(deaths_diff, 0), " (", deaths_diff_per, "%)")
oddeaths <- paste0(round(od_deaths_diff, 0), " (", od_deaths_diff_per, "%)")
effects_tbl <- cbind.data.frame(mod_names[-1], costs, deaths, oddeaths)
effects_tbl |>
flextable() |>
add_header_row(colwidths = c(1, 3),
values = c("", "Mean Change Compared to No Intervention")) |>
set_header_labels(values = list(
`mod_names[-1]` = "Interventions",
costs = "Discounted Net Present Costs \n in Millions (%)",
deaths = "Deaths (%)",
oddeaths = "Overdose-related Deaths (%)")) |>
theme_vanilla() |>
set_table_properties(layout = "autofit") |>
align_text_col(align = "center", header = TRUE, footer = TRUE)Mean Change Compared to No Intervention | |||
|---|---|---|---|
Interventions | Discounted Net Present Costs | Deaths (%) | Overdose-related Deaths (%) |
Naloxone | 102 (0.01%) | 396 (0.01%) | -6420 (-4.05%) |
Safer Supply | 4233 (0.53%) | -13956 (-0.31%) | -3138 (-1.98%) |
Prescription Guidelines | -26103 (-3.28%) | 15252 (0.34%) | -1764 (-1.11%) |
All Interventions | -24876 (-3.13%) | 14115 (0.31%) | -8543 (-5.39%) |
5.0.2 Scenario 2
Code
cost_diff_2 <- c(mod_nalx_2$total_net_present_cost - mod_basecase_2$total_net_present_cost,
mod_ss_2$total_net_present_cost - mod_basecase_2$total_net_present_cost,
mod_pres_guid_2$total_net_present_cost - mod_basecase_2$total_net_present_cost,
mod_all_int_2$total_net_present_cost - mod_basecase_2$total_net_present_cost)
cost_diff_per_2 <- c(round((cost_diff_2/c(mod_basecase_2$total_net_present_cost))*100,2))
deaths_eff_tbl_2 <- trace_tbl %>%
filter(scenario == "Scenario 2") %>%
filter(cycle_num == 180) %>%
select(-cycle_num) %>%
pivot_wider(names_from = mod, values_from = count) %>%
filter(state %in% c("BO_DEATH", "BO_OD_DEATH")) %>%
mutate(nalx_effect = Naloxone - `No Interventions`,
ss_effect = `Safer Supply` - `No Interventions`,
pg_effect = `Prescription Guidelines` - `No Interventions`,
all_effect = `All Interventions` - `No Interventions`,
nalx_effect_per = round(((Naloxone - `No Interventions`)/`No Interventions`)*100, 2),
ss_effect_per = round(((`Safer Supply` - `No Interventions`)/`No Interventions`)*100, 2),
pg_effect_per = round(((`Prescription Guidelines` - `No Interventions`)/`No Interventions`)*100, 2),
all_effect_per = round(((`All Interventions` - `No Interventions`)/`No Interventions`)*100, 2)) %>%
select(-c(Naloxone, `No Interventions`, `Safer Supply`, `Prescription Guidelines`, `All Interventions`))
od_deaths_diff_2 <- c(deaths_eff_tbl_2$nalx_effect[1], deaths_eff_tbl_2$ss_effect[1],
deaths_eff_tbl_2$pg_effect[1], deaths_eff_tbl_2$all_effect[1])
od_deaths_diff_per_2 <- c(deaths_eff_tbl_2$nalx_effect_per[1], deaths_eff_tbl_2$ss_effect_per[1],
deaths_eff_tbl_2$pg_effect_per[1], deaths_eff_tbl_2$all_effect_per[1])
deaths_diff_2 <- c(deaths_eff_tbl_2$nalx_effect[2], deaths_eff_tbl_2$ss_effect[2],
deaths_eff_tbl_2$pg_effect[2], deaths_eff_tbl_2$all_effect[2])
deaths_diff_per_2 <- c(deaths_eff_tbl_2$nalx_effect_per[2], deaths_eff_tbl_2$ss_effect_per[2],
deaths_eff_tbl_2$pg_effect_per[2], deaths_eff_tbl_2$all_effect_per[2])
costs_2 <- paste0(round(cost_diff_2/1000000, 0), " (", cost_diff_per_2, "%)")
deaths_2 <- paste0(round(deaths_diff_2, 0), " (", deaths_diff_per_2, "%)")
oddeaths_2 <- paste0(round(od_deaths_diff_2, 0), " (", od_deaths_diff_per_2, "%)")
effects_tbl_2 <- cbind.data.frame(mod_names[-1], costs_2, deaths_2, oddeaths_2)
effects_tbl_2 |>
flextable() |>
add_header_row(colwidths = c(1, 3),
values = c("", "Mean Change Compared to No Intervention")) |>
set_header_labels(values = list(
`mod_names[-1]` = "Interventions",
costs_2 = "Discounted Net Present Costs \n in Millions (%)",
deaths_2 = "Deaths (%)",
oddeaths_2 = "Overdose-related Deaths (%)")) |>
theme_vanilla() |>
set_table_properties(layout = "autofit") |>
align_text_col(align = "center", header = TRUE, footer = TRUE)Mean Change Compared to No Intervention | |||
|---|---|---|---|
Interventions | Discounted Net Present Costs | Deaths (%) | Overdose-related Deaths (%) |
Naloxone | 85 (0.01%) | 303 (0.01%) | -4861 (-3.9%) |
Safer Supply | 4047 (0.51%) | -12743 (-0.28%) | -2435 (-1.95%) |
Prescription Guidelines | -26098 (-3.29%) | 15308 (0.34%) | -1485 (-1.19%) |
All Interventions | -24882 (-3.13%) | 14105 (0.31%) | -6607 (-5.3%) |
Code
tbl_effect_plot1 <- cbind.data.frame(Intervention = mod_names[-1], cost_diff_per,
deaths_diff_per, od_deaths_diff_per) %>%
pivot_longer(2:4, names_to = "group", values_to = "percent_change") %>%
mutate(scenario = "Scenario 1") %>%
bind_rows(., cbind.data.frame(Intervention = mod_names[-1], cost_diff_per_2,
deaths_diff_per_2, od_deaths_diff_per_2) %>%
pivot_longer(2:4, names_to = "group", values_to = "percent_change") %>%
mutate(scenario = "Scenario 2")) %>%
mutate(group = ifelse(group == "cost_diff_per_2", "cost_diff_per",
ifelse(group == "deaths_diff_per_2", "deaths_diff_per",
ifelse(group == "od_deaths_diff_per_2", "od_deaths_diff_per", group))))5.0.2.1 Scenario 1
Code
p7 <- ggplot(data = tbl_effect_plot1 %>%
filter(scenario == "Scenario 1"), aes(x = group, y = percent_change))+
geom_point(aes(colour = Intervention, shape = scenario), size = 1.5) +
scale_color_brewer(palette = "Dark2") +
scale_shape_manual(values = c(15, 16)) +
geom_hline(yintercept = 0) +
ylab("Percent Change Compared to No Intervention (%)") +
xlab("") +
scale_x_discrete(name ="",
labels=c("Costs","Deaths","OD-releated Deaths"))+
theme_minimal()
ggplotly(p7)5.0.2.2 Scenario 2
Code
p8 <- ggplot(data = tbl_effect_plot1 %>%
filter(scenario == "Scenario 2"), aes(x = group, y = percent_change))+
geom_point(aes(colour = Intervention, shape = scenario), size = 1.5) +
scale_color_brewer(palette = "Dark2") +
scale_shape_manual(values = c(15, 16)) +
geom_hline(yintercept = 0) +
ylab("Percent Change Compared to No Intervention (%)") +
xlab("") +
scale_x_discrete(name ="",
labels=c("Costs","Deaths","OD-releated Deaths"))+
theme_minimal()
ggplotly(p8)5.0.2.3 Both
Code
ggplotly(ggplot(data = tbl_effect_plot1, aes(x = group, y = percent_change))+
geom_point(aes(colour = Intervention, shape = scenario), size = 1.5) +
scale_color_brewer(palette = "Dark2") +
# scale_shape_manual(values = c(15, 16, 17, 18)) +
scale_shape_manual(values = c(15, 16)) +
geom_hline(yintercept = 0) +
ylab("Percent Change Compared to No Intervention (%)") +
xlab("") +
scale_x_discrete(name ="",
labels=c("Costs","Deaths","OD-releated Deaths"))+
theme_bw() +
facet_wrap(~scenario))Code
tbl_effect_plot2 <- cbind.data.frame(Intervention = mod_names[-1], cost_diff_per,
od_deaths_diff_per) %>%
mutate(scenario = "Scenario 1") %>%
bind_rows(., cbind.data.frame(Intervention = mod_names[-1], cost_diff_per_2,
od_deaths_diff_per_2) %>%
rename(cost_diff_per = "cost_diff_per_2",
od_deaths_diff_per = "od_deaths_diff_per_2") %>%
mutate(scenario = "Scenario 2"))5.0.2.4 Scenario 1
Code
p9 <- ggplot(data = tbl_effect_plot2 %>% filter(scenario == "Scenario 1"),
aes(x = cost_diff_per, y = od_deaths_diff_per))+
geom_point(aes(color = Intervention, shape = scenario)) +
scale_color_brewer(palette = "Dark2") +
scale_shape_manual(values = c(15, 16)) +
geom_hline(yintercept = 0) +
geom_vline(xintercept = 0) +
xlab("Mean Change in Costs Compared to No Intervention (%)") +
ylab("Mean Change in Overdose-related deaths \n Compared to No Intervention (%)") +
theme_minimal()
ggplotly(p9)5.0.2.5 Scenario 2
Code
p10 <- ggplot(data = tbl_effect_plot2 %>% filter(scenario == "Scenario 2"),
aes(x = cost_diff_per, y = od_deaths_diff_per))+
geom_point(aes(color = Intervention, shape = scenario)) +
scale_color_brewer(palette = "Dark2") +
scale_shape_manual(values = c(15, 16)) +
geom_hline(yintercept = 0) +
geom_vline(xintercept = 0) +
xlab("Mean Change in Costs Compared to No Intervention (%)") +
ylab("Mean Change in Overdose-related deaths \n Compared to No Intervention (%)") +
theme_minimal()
ggplotly(p10)5.0.2.6 Both
Code
ggplotly(ggplot(data = tbl_effect_plot2, aes(x = cost_diff_per, y = od_deaths_diff_per))+
geom_point(aes(color = Intervention, shape = scenario)) +
scale_color_brewer(palette = "Dark2") +
scale_shape_manual(values = c(15, 16)) +
geom_hline(yintercept = 0) +
geom_vline(xintercept = 0) +
xlab("Mean Change in Costs Compared to No Intervention (%)") +
ylab("Mean Change in Overdose-related deaths \n Compared to No Intervention (%)") +
theme_bw() +
facet_wrap(~scenario))6 Secondary Outcomes
Code
v_yrs_prev <- c(2015:2029)
yrs_prev <- length(v_yrs_prev)
prev_fun <- function(mod, i, v_states){
prop_uncalib <- data.frame(matrix(NA,byrow = T, nrow = 12,
ncol = yrs_prev,
list(NULL,
paste0("prop_",
str_sub(v_yrs_prev, 3, 4)))))
tot_pop_uncalib_tbl <- data.frame(matrix(NA,byrow = T, nrow = 12,
ncol = yrs_prev,
list(NULL,
paste0("tot_pop_",
str_sub(v_yrs_prev, 3, 4)))))
for (j in 1:yrs_prev) {
yr <- (v_yrs_prev[1] - 1) + j
if (length(v_states) > 1)
{
prop_uncalib[, j] <- assign(
paste0("prop_", str_sub(yr, 3, 4)),
rowSums(mod$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr] + 1,
v_states]) /
rowSums(mod$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr] + 1,])
)
} else {
prop_uncalib[, j] <- assign(
paste0("prop_", str_sub(yr, 3, 4)),
mod$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr] + 1,
v_states]) /
rowSums(mod$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr] + 1,])
}
tot_pop_uncalib_tbl[, j] <- assign(
paste0("tot_pop_", str_sub(yr, 3, 4)),
rowSums(mod$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr] + 1,])
)
}
# final table for monthly prevalence over 15 years
prop_uncalib_tbl <- prop_uncalib %>%
pivot_longer(1:15, names_to = "grp", values_to = "value") %>%
arrange(grp) %>%
mutate(year = rep(v_yrs_prev,
each = 12),
year_mon = paste(year, month.abb[1:12], sep = "_")) %>%
bind_cols(., tot_pop_uncalib_tbl %>%
pivot_longer(1:15, names_to = "tot_pop",
values_to = "tot_pop_val") %>%
arrange(tot_pop) %>% select(-tot_pop)) %>%
mutate(intervention = mod_names[i])
# table with weighted yearly prevalence
prop_uncalib_wei_mean <- prop_uncalib_tbl %>%
group_by(year) %>%
summarise(value = weighted.mean(value, tot_pop_val)) %>%
ungroup() %>%
mutate(intervention = mod_names[i])
return (list(prop_uncalib_tbl = prop_uncalib_tbl,
prop_uncalib_wei_mean = prop_uncalib_wei_mean))
}6.1 Severe Brain injury
Code
noint_prev_sevbi_mon_tbl <- prev_fun(mod = mod_basecase, i = 1,
v_states = "BO_SEVERE_BI")$prop_uncalib_tbl
noint_prev_sevbi_yr_tbl <- prev_fun(mod = mod_basecase, i = 1,
v_states = "BO_SEVERE_BI")$prop_uncalib_wei_mean
nalx_prev_sevbi_mon_tbl <- prev_fun(mod = mod_nalx_1, i = 2,
v_states = "BO_SEVERE_BI")$prop_uncalib_tbl
nalx_prev_sevbi_yr_tbl <- prev_fun(mod = mod_nalx_1, i =2,
v_states = "BO_SEVERE_BI")$prop_uncalib_wei_mean
ss_prev_sevbi_mon_tbl <- prev_fun(mod = mod_ss_1, i = 3,
v_states = "BO_SEVERE_BI")$prop_uncalib_tbl
ss_prev_sevbi_yr_tbl <- prev_fun(mod = mod_ss_1, i =3,
v_states = "BO_SEVERE_BI")$prop_uncalib_wei_mean
pg_prev_sevbi_mon_tbl <- prev_fun(mod = mod_pres_guid_1, i = 4,
v_states = "BO_SEVERE_BI")$prop_uncalib_tbl
pg_prev_sevbi_yr_tbl <- prev_fun(mod = mod_pres_guid_1, i =4,
v_states = "BO_SEVERE_BI")$prop_uncalib_wei_mean
allint_prev_sevbi_mon_tbl <- prev_fun(mod = mod_all_int_1, i = 5,
v_states = "BO_SEVERE_BI")$prop_uncalib_tbl
allint_prev_sevbi_yr_tbl <- prev_fun(mod = mod_all_int_1, i =5,
v_states = "BO_SEVERE_BI")$prop_uncalib_wei_mean
prev_sevbi_mon_tbl <- noint_prev_sevbi_mon_tbl %>%
bind_rows(., nalx_prev_sevbi_mon_tbl, ss_prev_sevbi_mon_tbl,
pg_prev_sevbi_mon_tbl, allint_prev_sevbi_mon_tbl)
prev_sevbi_yr_tbl <- noint_prev_sevbi_yr_tbl %>%
bind_rows(., nalx_prev_sevbi_yr_tbl, ss_prev_sevbi_yr_tbl,
pg_prev_sevbi_yr_tbl, allint_prev_sevbi_yr_tbl)
prev_sevbi_mon_tbl$year_mon <- factor(prev_sevbi_mon_tbl$year_mon,
levels = rep(
paste(
rep(v_yrs_prev,
each = 12),
month.abb[1:12],
sep = "_"
), 5),
labels = rep(
paste(
rep(v_yrs_prev,
each = 12),
month.abb[1:12],
sep = "_"), 5))
ggplotly(ggplot(trace_tbl %>%
filter(state %in% "BO_SEVERE_BI") %>%
filter(scenario == "Scenario 1") %>%
left_join(., year_mon_cycle_tbl, by = c("cycle_num" = "cycle")),
aes(x = cycle_num, y = count, colour = mod, key = year)) +
geom_line() +
ylab("Count") + xlab("Cycle Number") +
scale_color_manual(values = colors_pal, name = "Intervention"))Code
ggplotly(ggplot() +
geom_point(data = prev_sevbi_yr_tbl,
aes(x = factor(year), y = value,
color = intervention, fill = intervention)) +
xlab("Year") + ylab("Weighted Yearly Prevalence of severe brain injury") +
scale_fill_discrete(name = "Intervention")+
scale_color_discrete(name = "Intervention")+
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)))6.2 Prevalence of moderate Brain injury
Code
noint_prev_modbi_mon_tbl <- prev_fun(mod = mod_basecase, i = 1,
v_states = "BO_MOD_BI")$prop_uncalib_tbl
noint_prev_modbi_yr_tbl <- prev_fun(mod = mod_basecase, i = 1,
v_states = "BO_MOD_BI")$prop_uncalib_wei_mean
nalx_prev_modbi_mon_tbl <- prev_fun(mod = mod_nalx_1, i = 2,
v_states = "BO_MOD_BI")$prop_uncalib_tbl
nalx_prev_modbi_yr_tbl <- prev_fun(mod = mod_nalx_1, i =2,
v_states = "BO_MOD_BI")$prop_uncalib_wei_mean
ss_prev_modbi_mon_tbl <- prev_fun(mod = mod_ss_1, i = 3,
v_states = "BO_MOD_BI")$prop_uncalib_tbl
ss_prev_modbi_yr_tbl <- prev_fun(mod = mod_ss_1, i =3,
v_states = "BO_MOD_BI")$prop_uncalib_wei_mean
pg_prev_modbi_mon_tbl <- prev_fun(mod = mod_pres_guid_1, i = 4,
v_states = "BO_MOD_BI")$prop_uncalib_tbl
pg_prev_modbi_yr_tbl <- prev_fun(mod = mod_pres_guid_1, i =4,
v_states = "BO_MOD_BI")$prop_uncalib_wei_mean
allint_prev_modbi_mon_tbl <- prev_fun(mod = mod_all_int_1, i = 5,
v_states = "BO_MOD_BI")$prop_uncalib_tbl
allint_prev_modbi_yr_tbl <- prev_fun(mod = mod_all_int_1, i =5,
v_states = "BO_MOD_BI")$prop_uncalib_wei_mean
prev_modbi_mon_tbl <- noint_prev_modbi_mon_tbl %>%
bind_rows(., nalx_prev_modbi_mon_tbl, ss_prev_modbi_mon_tbl,
pg_prev_modbi_mon_tbl, allint_prev_modbi_mon_tbl)
prev_modbi_yr_tbl <- noint_prev_modbi_yr_tbl %>%
bind_rows(., nalx_prev_modbi_yr_tbl, ss_prev_modbi_yr_tbl,
pg_prev_modbi_yr_tbl, allint_prev_modbi_yr_tbl)
prev_modbi_mon_tbl$year_mon <- factor(prev_modbi_mon_tbl$year_mon,
levels = rep(
paste(
rep(v_yrs_prev,
each = 12),
month.abb[1:12],
sep = "_"
), 5),
labels = rep(
paste(
rep(v_yrs_prev,
each = 12),
month.abb[1:12],
sep = "_"), 5))
ggplotly(ggplot(trace_tbl %>%
filter(state %in% "BO_MOD_BI") %>%
filter(scenario == "Scenario 1") %>%
group_by(cycle_num, mod) %>%
summarize(count = sum(count)) %>%
ungroup() %>%
left_join(., year_mon_cycle_tbl, by = c("cycle_num" = "cycle")),
aes(x = cycle_num, y = count, colour = mod, key = year)) +
geom_line() +
ylab("Count") + xlab("Cycle Number") +
scale_color_manual(values = colors_pal, name = "Intervention"))Code
ggplotly(ggplot() +
geom_point(data = prev_modbi_yr_tbl,
aes(x = factor(year), y = value,
color = intervention, fill = intervention)) +
xlab("Year") + ylab("Weighted Yearly Prevalence of moderate brain injury") +
scale_fill_discrete(name = "Intervention")+
scale_color_discrete(name = "Intervention")+
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)))6.3 Prevalence of substance use treatment
Code
v_states_prev_tx <- c("BS_DETOX", "BS_OAT_INI", "BS_OAT_MAINT",
"BS_OAT_SS", "BS_SS", "BR_OAT_INI", "BR_OAT_MAINT",
"BR_OAT_SS", "BR_SS")
noint_prev_tx_mon_tbl <- prev_fun(mod = mod_basecase, i = 1,
v_states = v_states_prev_tx)$prop_uncalib_tbl
noint_prev_tx_yr_tbl <- prev_fun(mod = mod_basecase, i = 1,
v_states = v_states_prev_tx)$prop_uncalib_wei_mean
nalx_prev_tx_mon_tbl <- prev_fun(mod = mod_nalx_1, i = 2,
v_states = v_states_prev_tx)$prop_uncalib_tbl
nalx_prev_tx_yr_tbl <- prev_fun(mod = mod_nalx_1, i =2,
v_states = v_states_prev_tx)$prop_uncalib_wei_mean
ss_prev_tx_mon_tbl <- prev_fun(mod = mod_ss_1, i = 3,
v_states = v_states_prev_tx)$prop_uncalib_tbl
ss_prev_tx_yr_tbl <- prev_fun(mod = mod_ss_1, i =3,
v_states = v_states_prev_tx)$prop_uncalib_wei_mean
pg_prev_tx_mon_tbl <- prev_fun(mod = mod_pres_guid_1, i = 4,
v_states = v_states_prev_tx)$prop_uncalib_tbl
pg_prev_tx_yr_tbl <- prev_fun(mod = mod_pres_guid_1, i =4,
v_states = v_states_prev_tx)$prop_uncalib_wei_mean
allint_prev_tx_mon_tbl <- prev_fun(mod = mod_all_int_1, i = 5,
v_states = v_states_prev_tx)$prop_uncalib_tbl
allint_prev_tx_yr_tbl <- prev_fun(mod = mod_all_int_1, i =5,
v_states = v_states_prev_tx)$prop_uncalib_wei_mean
prev_tx_mon_tbl <- noint_prev_tx_mon_tbl %>%
bind_rows(., nalx_prev_tx_mon_tbl, ss_prev_tx_mon_tbl,
pg_prev_tx_mon_tbl, allint_prev_tx_mon_tbl)
prev_tx_yr_tbl <- noint_prev_tx_yr_tbl %>%
bind_rows(., nalx_prev_tx_yr_tbl, ss_prev_tx_yr_tbl,
pg_prev_tx_yr_tbl, allint_prev_tx_yr_tbl)
prev_tx_mon_tbl$year_mon <- factor(prev_tx_mon_tbl$year_mon,
levels = rep(
paste(
rep(v_yrs_prev,
each = 12),
month.abb[1:12],
sep = "_"
), 5),
labels = rep(
paste(
rep(v_yrs_prev,
each = 12),
month.abb[1:12],
sep = "_"), 5))
ggplotly(ggplot(trace_tbl %>%
filter(state %in% v_states_prev_tx) %>%
filter(scenario == "Scenario 1") %>%
group_by(cycle_num, mod) %>%
summarize(count = sum(count)) %>%
ungroup() %>%
left_join(., year_mon_cycle_tbl, by = c("cycle_num" = "cycle")),
aes(x = cycle_num, y = count, colour = mod, key = year)) +
geom_line() +
ylab("Count") + xlab("Cycle Number") +
scale_color_manual(values = colors_pal, name = "Intervention"))Code
ggplotly(ggplot() +
geom_point(data = prev_tx_yr_tbl,
aes(x = factor(year), y = value,
color = intervention, fill = intervention)) +
xlab("Year") + ylab("Weighted Yearly Prevalence of substance use treatment") +
scale_fill_discrete(name = "Intervention")+
scale_color_discrete(name = "Intervention")+
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)))6.4 Prevalence of prescription opioid use
Code
prop_opioids_rx_target_tbl <- calib_target_tbl %>%
filter(group == "prev_on_opioidsrx") %>%
select(year, target) %>%
rename(value = target) %>%
mutate(year_mon = paste(year, month.abb[7], sep = "_"))
v_states_prev_rx_opd <- c("BPO_ACUTE", "BPO_CHRONIC", "BPO_CANCER",
"BPO_OTHER", "BPO_PALLIATIVE")
noint_prev_rx_opd_mon_tbl <- prev_fun(mod = mod_basecase, i = 1,
v_states = v_states_prev_rx_opd)$prop_uncalib_tbl
noint_prev_rx_opd_yr_tbl <- prev_fun(mod = mod_basecase, i = 1,
v_states = v_states_prev_rx_opd)$prop_uncalib_wei_mean
nalx_prev_rx_opd_mon_tbl <- prev_fun(mod = mod_nalx_1, i = 2,
v_states = v_states_prev_rx_opd)$prop_uncalib_tbl
nalx_prev_rx_opd_yr_tbl <- prev_fun(mod = mod_nalx_1, i =2,
v_states = v_states_prev_rx_opd)$prop_uncalib_wei_mean
ss_prev_rx_opd_mon_tbl <- prev_fun(mod = mod_ss_1, i = 3,
v_states = v_states_prev_rx_opd)$prop_uncalib_tbl
ss_prev_rx_opd_yr_tbl <- prev_fun(mod = mod_ss_1, i =3,
v_states = v_states_prev_rx_opd)$prop_uncalib_wei_mean
pg_prev_rx_opd_mon_tbl <- prev_fun(mod = mod_pres_guid_1, i = 4,
v_states = v_states_prev_rx_opd)$prop_uncalib_tbl
pg_prev_rx_opd_yr_tbl <- prev_fun(mod = mod_pres_guid_1, i =4,
v_states = v_states_prev_rx_opd)$prop_uncalib_wei_mean
allint_prev_rx_opd_mon_tbl <- prev_fun(mod = mod_all_int_1, i = 5,
v_states = v_states_prev_rx_opd)$prop_uncalib_tbl
allint_prev_rx_opd_yr_tbl <- prev_fun(mod = mod_all_int_1, i =5,
v_states = v_states_prev_rx_opd)$prop_uncalib_wei_mean
prev_rx_opd_mon_tbl <- noint_prev_rx_opd_mon_tbl %>%
bind_rows(., nalx_prev_rx_opd_mon_tbl, ss_prev_rx_opd_mon_tbl,
pg_prev_rx_opd_mon_tbl, allint_prev_rx_opd_mon_tbl) %>%
bind_rows(., prop_opioids_rx_target_tbl %>%
select(year, value, year_mon) %>%
mutate(intervention = "Target"))
prev_rx_opd_yr_tbl <- noint_prev_rx_opd_yr_tbl %>%
bind_rows(., nalx_prev_rx_opd_yr_tbl, ss_prev_rx_opd_yr_tbl,
pg_prev_rx_opd_yr_tbl, allint_prev_rx_opd_yr_tbl) %>%
bind_rows(., prop_opioids_rx_target_tbl %>%
select(year, value) %>%
mutate(intervention = "Target"))
prev_rx_opd_mon_tbl$year_mon <- factor(prev_rx_opd_mon_tbl$year_mon,
levels = rep(
paste(
rep(v_yrs_prev,
each = 12),
month.abb[1:12],
sep = "_"
), 5),
labels = rep(
paste(
rep(v_yrs_prev,
each = 12),
month.abb[1:12],
sep = "_"), 5))
prop_opioids_rx_target_tbl$year_mon <- factor(prop_opioids_rx_target_tbl$year_mon,
levels = paste(rep(v_yrs_prev,
each = 12),
month.abb[1:12],
sep = "_"),
labels = paste(rep(v_yrs_prev,
each = 12),
month.abb[1:12],
sep = "_"))
ggplotly(ggplot(trace_tbl %>%
filter(state %in% c("BPO_ACUTE", "BPO_CHRONIC", "BPO_CANCER",
"BPO_OTHER", "BPO_PALLIATIVE")) %>%
filter(scenario == "Scenario 1") %>%
group_by(cycle_num, mod) %>%
summarize(count = sum(count)) %>%
ungroup() %>%
left_join(., year_mon_cycle_tbl, by = c("cycle_num" = "cycle")),
aes(x = cycle_num, y = count, colour = mod, key = year)) +
geom_line() +
ylab("Count") + xlab("Cycle Number") +
scale_color_manual(values = colors_pal, name = "Intervention"))Code
ggplotly(ggplot() +
geom_point(data = prev_rx_opd_yr_tbl,
aes(x = factor(year), y = value,
color = intervention, fill = intervention)) +
xlab("Year") + ylab("Weighted Yearly Prevalence of prescription opioid use") +
scale_fill_discrete(name = "Intervention")+
scale_color_discrete(name = "Intervention")+
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)))6.5 Prevalence of illicit use
Code
v_states_illicit <- c("BI_ILLICIT", "BR_ILLICIT")
noint_prev_illicituse_mon_tbl <- prev_fun(mod = mod_basecase, i = 1,
v_states = v_states_illicit)$prop_uncalib_tbl
noint_prev_illicituse_yr_tbl <- prev_fun(mod = mod_basecase, i = 1,
v_states = v_states_illicit)$prop_uncalib_wei_mean
nalx_prev_illicituse_mon_tbl <- prev_fun(mod = mod_nalx_1, i = 2,
v_states = v_states_illicit)$prop_uncalib_tbl
nalx_prev_illicituse_yr_tbl <- prev_fun(mod = mod_nalx_1, i =2,
v_states = v_states_illicit)$prop_uncalib_wei_mean
ss_prev_illicituse_mon_tbl <- prev_fun(mod = mod_ss_1, i = 3,
v_states = v_states_illicit)$prop_uncalib_tbl
ss_prev_illicituse_yr_tbl <- prev_fun(mod = mod_ss_1, i =3,
v_states = v_states_illicit)$prop_uncalib_wei_mean
pg_prev_illicituse_mon_tbl <- prev_fun(mod = mod_pres_guid_1, i = 4,
v_states = v_states_illicit)$prop_uncalib_tbl
pg_prev_illicituse_yr_tbl <- prev_fun(mod = mod_pres_guid_1, i =4,
v_states = v_states_illicit)$prop_uncalib_wei_mean
allint_prev_illicituse_mon_tbl <- prev_fun(mod = mod_all_int_1, i = 5,
v_states = v_states_illicit)$prop_uncalib_tbl
allint_prev_illicituse_yr_tbl <- prev_fun(mod = mod_all_int_1, i =5,
v_states = v_states_illicit)$prop_uncalib_wei_mean
prev_illicituse_mon_tbl <- noint_prev_illicituse_mon_tbl %>%
bind_rows(., nalx_prev_illicituse_mon_tbl, ss_prev_illicituse_mon_tbl,
pg_prev_illicituse_mon_tbl, allint_prev_illicituse_mon_tbl)
prev_illicituse_yr_tbl <- noint_prev_illicituse_yr_tbl %>%
bind_rows(., nalx_prev_illicituse_yr_tbl, ss_prev_illicituse_yr_tbl,
pg_prev_illicituse_yr_tbl, allint_prev_illicituse_yr_tbl)
prev_illicituse_mon_tbl$year_mon <- factor(prev_illicituse_mon_tbl$year_mon,
levels = rep(
paste(
rep(v_yrs_prev,
each = 12),
month.abb[1:12],
sep = "_"
), 5),
labels = rep(
paste(
rep(v_yrs_prev,
each = 12),
month.abb[1:12],
sep = "_"), 5))
ggplotly(ggplot(trace_tbl %>%
filter(state %in% v_states_illicit) %>%
filter(scenario == "Scenario 1") %>%
group_by(cycle_num, mod) %>%
summarize(count = sum(count)) %>%
ungroup() %>%
left_join(., year_mon_cycle_tbl, by = c("cycle_num" = "cycle")),
aes(x = cycle_num, y = count, colour = mod, key = year)) +
geom_line() +
ylab("Count") + xlab("Cycle Number") +
scale_color_manual(values = colors_pal, name = "Intervention"))Code
ggplotly(ggplot() +
geom_point(data = prev_illicituse_yr_tbl,
aes(x = factor(year), y = value,
color = intervention, fill = intervention)) +
xlab("Year") + ylab("Weighted Yearly Prevalence of illicit opioid use") +
scale_fill_discrete(name = "Intervention")+
scale_color_discrete(name = "Intervention")+
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)))6.6 Prevalence of misuse
Code
noint_prev_rxmisuse_mon_tbl <- prev_fun(mod = mod_basecase, i = 1,
v_states = "BPO_MISUSE")$prop_uncalib_tbl
noint_prev_rxmisuse_yr_tbl <- prev_fun(mod = mod_basecase, i = 1,
v_states = "BPO_MISUSE")$prop_uncalib_wei_mean
nalx_prev_rxmisuse_mon_tbl <- prev_fun(mod = mod_nalx_1, i = 2,
v_states = "BPO_MISUSE")$prop_uncalib_tbl
nalx_prev_rxmisuse_yr_tbl <- prev_fun(mod = mod_nalx_1, i =2,
v_states = "BPO_MISUSE")$prop_uncalib_wei_mean
ss_prev_rxmisuse_mon_tbl <- prev_fun(mod = mod_ss_1, i = 3,
v_states = "BPO_MISUSE")$prop_uncalib_tbl
ss_prev_rxmisuse_yr_tbl <- prev_fun(mod = mod_ss_1, i =3,
v_states = "BPO_MISUSE")$prop_uncalib_wei_mean
pg_prev_rxmisuse_mon_tbl <- prev_fun(mod = mod_pres_guid_1, i = 4,
v_states = "BPO_MISUSE")$prop_uncalib_tbl
pg_prev_rxmisuse_yr_tbl <- prev_fun(mod = mod_pres_guid_1, i =4,
v_states = "BPO_MISUSE")$prop_uncalib_wei_mean
allint_prev_rxmisuse_mon_tbl <- prev_fun(mod = mod_all_int_1, i = 5,
v_states = "BPO_MISUSE")$prop_uncalib_tbl
allint_prev_rxmisuse_yr_tbl <- prev_fun(mod = mod_all_int_1, i =5,
v_states = "BPO_MISUSE")$prop_uncalib_wei_mean
prev_rxmisuse_mon_tbl <- noint_prev_rxmisuse_mon_tbl %>%
bind_rows(., nalx_prev_rxmisuse_mon_tbl, ss_prev_rxmisuse_mon_tbl,
pg_prev_rxmisuse_mon_tbl, allint_prev_rxmisuse_mon_tbl)
prev_rxmisuse_yr_tbl <- noint_prev_rxmisuse_yr_tbl %>%
bind_rows(., nalx_prev_rxmisuse_yr_tbl, ss_prev_rxmisuse_yr_tbl,
pg_prev_rxmisuse_yr_tbl, allint_prev_rxmisuse_yr_tbl)
prev_rxmisuse_mon_tbl$year_mon <- factor(prev_rxmisuse_mon_tbl$year_mon,
levels = rep(
paste(
rep(v_yrs_prev,
each = 12),
month.abb[1:12],
sep = "_"
), 5),
labels = rep(
paste(
rep(v_yrs_prev,
each = 12),
month.abb[1:12],
sep = "_"), 5))
ggplotly(ggplot(trace_tbl %>%
filter(state %in% "BPO_MISUSE") %>%
filter(scenario == "Scenario 1") %>%
group_by(cycle_num, mod) %>%
summarize(count = sum(count)) %>%
ungroup() %>%
left_join(., year_mon_cycle_tbl, by = c("cycle_num" = "cycle")),
aes(x = cycle_num, y = count, colour = mod, key = year)) +
geom_line() +
ylab("Count") + xlab("Cycle Number") +
scale_color_manual(values = colors_pal, name = "Intervention"))Code
ggplotly(ggplot() +
geom_point(data = prev_rxmisuse_yr_tbl,
aes(x = factor(year), y = value,
color = intervention, fill = intervention)) +
xlab("Year") + ylab("Weighted Yearly Prevalence of rx opioid misuse") +
scale_fill_discrete(name = "Intervention")+
scale_color_discrete(name = "Intervention")+
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)))6.7 Prevalence of overdose
Code
v_states_od <- c("BO_OD_RX", "BO_OD_ILLICIT", "BR_OD_ILLICIT")
noint_prev_od_mon_tbl <- prev_fun(mod = mod_basecase, i = 1,
v_states = v_states_od)$prop_uncalib_tbl
noint_prev_od_yr_tbl <- prev_fun(mod = mod_basecase, i = 1,
v_states = v_states_od)$prop_uncalib_wei_mean
nalx_prev_od_mon_tbl <- prev_fun(mod = mod_nalx_1, i = 2,
v_states = v_states_od)$prop_uncalib_tbl
nalx_prev_od_yr_tbl <- prev_fun(mod = mod_nalx_1, i =2,
v_states = v_states_od)$prop_uncalib_wei_mean
ss_prev_od_mon_tbl <- prev_fun(mod = mod_ss_1, i = 3,
v_states = v_states_od)$prop_uncalib_tbl
ss_prev_od_yr_tbl <- prev_fun(mod = mod_ss_1, i =3,
v_states = v_states_od)$prop_uncalib_wei_mean
pg_prev_od_mon_tbl <- prev_fun(mod = mod_pres_guid_1, i = 4,
v_states = v_states_od)$prop_uncalib_tbl
pg_prev_od_yr_tbl <- prev_fun(mod = mod_pres_guid_1, i =4,
v_states = v_states_od)$prop_uncalib_wei_mean
allint_prev_od_mon_tbl <- prev_fun(mod = mod_all_int_1, i = 5,
v_states = v_states_od)$prop_uncalib_tbl
allint_prev_od_yr_tbl <- prev_fun(mod = mod_all_int_1, i =5,
v_states = v_states_od)$prop_uncalib_wei_mean
prev_od_mon_tbl <- noint_prev_od_mon_tbl %>%
bind_rows(., nalx_prev_od_mon_tbl, ss_prev_od_mon_tbl,
pg_prev_od_mon_tbl, allint_prev_od_mon_tbl)
prev_od_yr_tbl <- noint_prev_od_yr_tbl %>%
bind_rows(., nalx_prev_od_yr_tbl, ss_prev_od_yr_tbl,
pg_prev_od_yr_tbl, allint_prev_od_yr_tbl)
prev_od_mon_tbl$year_mon <- factor(prev_od_mon_tbl$year_mon,
levels = rep(
paste(
rep(v_yrs_prev,
each = 12),
month.abb[1:12],
sep = "_"
), 5),
labels = rep(
paste(
rep(v_yrs_prev,
each = 12),
month.abb[1:12],
sep = "_"), 5))
ggplotly(ggplot(trace_tbl %>%
filter(state %in% v_states_od) %>%
filter(scenario == "Scenario 1") %>%
group_by(cycle_num, mod) %>%
summarize(count = sum(count)) %>%
ungroup() %>%
left_join(., year_mon_cycle_tbl, by = c("cycle_num" = "cycle")),
aes(x = cycle_num, y = count, colour = mod, key = year)) +
geom_line() +
ylab("Count") + xlab("Cycle Number") +
scale_color_manual(values = colors_pal, name = "Intervention"))Code
ggplotly(ggplot() +
geom_point(data = prev_od_yr_tbl,
aes(x = factor(year), y = value,
color = intervention, fill = intervention)) +
xlab("Year") + ylab("Weighted Yearly Prevalence of opioid-related overdose") +
scale_fill_discrete(name = "Intervention")+
scale_color_discrete(name = "Intervention")+
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)))6.7.1 prevalence of illicit overdose
Code
v_states_illicitod <- c("BO_OD_ILLICIT", "BR_OD_ILLICIT")
noint_prev_illicitod_mon_tbl <- prev_fun(mod = mod_basecase, i = 1,
v_states = v_states_illicitod)$prop_uncalib_tbl
noint_prev_illicitod_yr_tbl <- prev_fun(mod = mod_basecase, i = 1,
v_states = v_states_illicitod)$prop_uncalib_wei_mean
nalx_prev_illicitod_mon_tbl <- prev_fun(mod = mod_nalx_1, i = 2,
v_states = v_states_illicitod)$prop_uncalib_tbl
nalx_prev_illicitod_yr_tbl <- prev_fun(mod = mod_nalx_1, i =2,
v_states = v_states_illicitod)$prop_uncalib_wei_mean
ss_prev_illicitod_mon_tbl <- prev_fun(mod = mod_ss_1, i = 3,
v_states = v_states_illicitod)$prop_uncalib_tbl
ss_prev_illicitod_yr_tbl <- prev_fun(mod = mod_ss_1, i =3,
v_states = v_states_illicitod)$prop_uncalib_wei_mean
pg_prev_illicitod_mon_tbl <- prev_fun(mod = mod_pres_guid_1, i = 4,
v_states = v_states_illicitod)$prop_uncalib_tbl
pg_prev_illicitod_yr_tbl <- prev_fun(mod = mod_pres_guid_1, i =4,
v_states = v_states_illicitod)$prop_uncalib_wei_mean
allint_prev_illicitod_mon_tbl <- prev_fun(mod = mod_all_int_1, i = 5,
v_states = v_states_illicitod)$prop_uncalib_tbl
allint_prev_illicitod_yr_tbl <- prev_fun(mod = mod_all_int_1, i =5,
v_states = v_states_illicitod)$prop_uncalib_wei_mean
prev_illicitod_mon_tbl <- noint_prev_illicitod_mon_tbl %>%
bind_rows(., nalx_prev_illicitod_mon_tbl, ss_prev_illicitod_mon_tbl,
pg_prev_illicitod_mon_tbl, allint_prev_illicitod_mon_tbl)
prev_illicitod_yr_tbl <- noint_prev_illicitod_yr_tbl %>%
bind_rows(., nalx_prev_illicitod_yr_tbl, ss_prev_illicitod_yr_tbl,
pg_prev_illicitod_yr_tbl, allint_prev_illicitod_yr_tbl)
prev_illicitod_mon_tbl$year_mon <- factor(prev_illicitod_mon_tbl$year_mon,
levels = rep(
paste(
rep(v_yrs_prev,
each = 12),
month.abb[1:12],
sep = "_"
), 5),
labels = rep(
paste(
rep(v_yrs_prev,
each = 12),
month.abb[1:12],
sep = "_"), 5))
ggplotly(ggplot(trace_tbl %>%
filter(state %in% v_states_illicitod) %>%
filter(scenario == "Scenario 1") %>%
group_by(cycle_num, mod) %>%
summarize(count = sum(count)) %>%
ungroup() %>%
left_join(., year_mon_cycle_tbl, by = c("cycle_num" = "cycle")),
aes(x = cycle_num, y = count, colour = mod, key = year)) +
geom_line() +
ylab("Count") + xlab("Cycle Number") +
scale_color_manual(values = colors_pal, name = "Intervention"))Code
ggplotly(ggplot() +
geom_point(data = prev_illicitod_yr_tbl,
aes(x = factor(year), y = value,
color = intervention, fill = intervention)) +
xlab("Year") + ylab("Weighted Yearly Prevalence of illicit opioid-related overdose") +
scale_fill_discrete(name = "Intervention")+
scale_color_discrete(name = "Intervention")+
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)))6.7.2 prevalence of rx overdose
Code
noint_prev_rxod_mon_tbl <- prev_fun(mod = mod_basecase, i = 1,
v_states = "BO_OD_RX")$prop_uncalib_tbl
noint_prev_rxod_yr_tbl <- prev_fun(mod = mod_basecase, i = 1,
v_states = "BO_OD_RX")$prop_uncalib_wei_mean
nalx_prev_rxod_mon_tbl <- prev_fun(mod = mod_nalx_1, i = 2,
v_states = "BO_OD_RX")$prop_uncalib_tbl
nalx_prev_rxod_yr_tbl <- prev_fun(mod = mod_nalx_1, i =2,
v_states = "BO_OD_RX")$prop_uncalib_wei_mean
ss_prev_rxod_mon_tbl <- prev_fun(mod = mod_ss_1, i = 3,
v_states = "BO_OD_RX")$prop_uncalib_tbl
ss_prev_rxod_yr_tbl <- prev_fun(mod = mod_ss_1, i =3,
v_states = "BO_OD_RX")$prop_uncalib_wei_mean
pg_prev_rxod_mon_tbl <- prev_fun(mod = mod_pres_guid_1, i = 4,
v_states = "BO_OD_RX")$prop_uncalib_tbl
pg_prev_rxod_yr_tbl <- prev_fun(mod = mod_pres_guid_1, i =4,
v_states = "BO_OD_RX")$prop_uncalib_wei_mean
allint_prev_rxod_mon_tbl <- prev_fun(mod = mod_all_int_1, i = 5,
v_states = "BO_OD_RX")$prop_uncalib_tbl
allint_prev_rxod_yr_tbl <- prev_fun(mod = mod_all_int_1, i =5,
v_states = "BO_OD_RX")$prop_uncalib_wei_mean
prev_rxod_mon_tbl <- noint_prev_rxod_mon_tbl %>%
bind_rows(., nalx_prev_rxod_mon_tbl, ss_prev_rxod_mon_tbl,
pg_prev_rxod_mon_tbl, allint_prev_rxod_mon_tbl)
prev_rxod_yr_tbl <- noint_prev_rxod_yr_tbl %>%
bind_rows(., nalx_prev_rxod_yr_tbl, ss_prev_rxod_yr_tbl,
pg_prev_rxod_yr_tbl, allint_prev_rxod_yr_tbl)
prev_rxod_mon_tbl$year_mon <- factor(prev_rxod_mon_tbl$year_mon,
levels = rep(
paste(
rep(v_yrs_prev,
each = 12),
month.abb[1:12],
sep = "_"
), 5),
labels = rep(
paste(
rep(v_yrs_prev,
each = 12),
month.abb[1:12],
sep = "_"), 5))
ggplotly(ggplot(trace_tbl %>%
filter(state %in% "BO_OD_RX") %>%
filter(scenario == "Scenario 1") %>%
group_by(cycle_num, mod) %>%
summarize(count = sum(count)) %>%
ungroup() %>%
left_join(., year_mon_cycle_tbl, by = c("cycle_num" = "cycle")),
aes(x = cycle_num, y = count, colour = mod, key = year)) +
geom_line() +
ylab("Count") + xlab("Cycle Number") +
scale_color_manual(values = colors_pal, name = "Intervention"))Code
ggplotly(ggplot() +
geom_point(data = prev_rxod_yr_tbl,
aes(x = factor(year), y = value,
color = intervention, fill = intervention)) +
xlab("Year") + ylab("Weighted Yearly Prevalence of rx opioid-related overdose") +
scale_fill_discrete(name = "Intervention")+
scale_color_discrete(name = "Intervention")+
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)))